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We propose to use the resonant enhancement of the parametric instability in an optomechanical system of 
two optical modes coupled to a mechanical oscillator to prepare mechanical states with sub-Poissonian phonon 
statistics. Strong single photon coupling is not required. The requirements regarding sideband resolution, 
circulating cavity power and environmental temperature are in reach with state of the art parameters of optome¬ 
chanical crystals. Phonon antibunching can be verfied in a Hanburry-Brown-Twiss measurement on the output 
field of the optomechanical cavity. 


Introduction Optomechanical experiments, where light 
resonators are coupled to mechanical oscillators (H [2|, are 
achieving increasingly good control of macroscopic objects 
on the quantum level: Milestones such as cooling the motion 
of these oscillators to their quantum ground state EUD, co¬ 
herent transfer of information between light and mirror |[5][6l, 
observation of radiation pressure shot noise on the oscillator 
(TllHl, as well as entanglement between the light field and the 
mirrors O have been achieved in recent years. 

The phonon analogue of a laser, which is realized using 
the optical cavity as the gain medium to excite coherent os¬ 
cillations of the mechanical oscillator has been demonstrated 
in (ToHISl, and its phonon statistics has been mapped out 
via a Hanburry-Brown-Twiss measurement on the sideband- 
photons emitted from the optomechanical cavity CD. The- 
oretical work suggests that it is possible to prepare a state 
with quantum signatures in the phonon statistics such as 
phonon antibunching and even negative Wigner density CS- 
[24ll . However, the requirements on system parameters to 
see phonon antibunching scale unfavorably, so that sub- 
Poissonian phonon statistics has eluded experimental obser¬ 
vation. 

In this article we propose to make use of the enhanced op¬ 
tomechanical nonlinearity 125^271 of a setup with two op¬ 
tical modes to overcome this difficulty and prepare phonon 
laser states featuring antibunching in steady state with state 
of the art optomechanical crystals. The enhanced nonlinear¬ 
ity has been discussed in the context of detectors for phonons 
or photons [261, quantum memory (281, and to improve (271 
the parameters of mechanically induced photon antibunching 
(29j |30l . In the context of the phonon laser transition the en¬ 
hanced optomechanical instability with two optical modes has 
been anticipated as a possible complication for gravitational 
wave detectors ED , and has been studied experimentally El- 
[BiEa and theoretically (33H361 in the classical regime. Here 
we show for the first time that one can detect quantum sig¬ 
natures in the phonon lasing of such a three-mode system. 
In particular phonon antibunching and, with more demand¬ 
ing system requirements, negative mechanical Wigner density 
can be prepared in steady state. 

Terminology for phonon statistics Denoting the phonon 
number h = c^c, its statistics is characterized by the Fano fac¬ 



Figure 1. left: Two optical modes a and b are coupled to a mechanical 
mode c. The b-modc is resonantly driven by a laser of strength E and 
frequency cjl - cot- The a-mode is detuned from b by A = cOm, the 
mechanical resonance frequency, as depicted in the plot on the right. 
The nonlinear interaction of the three modes a, b, and c gives rise 
to optomechanical limit cycles with strongly sub-Poissonian phonon 
number statistics. A third optical mode d can be used to reduce the 
effective temperature of the mechanical oscillator’s bath and to read 
out the phonon statistics in a Hanburry-Brown-Twiss measurement. 


tor F = {Ah^)l{h), and the second order coherence function 
g^^\t) at time t = 0 

^(2)(Q) ^ = l+(F- (1) 

{ny 

which gives information on the temporal correlations of the 
phonons. ( g^^\0) > 1 and < 1 corresponding to 

bunching and anti-bunching respectively E3-) The Fano 
factor F can be inferred from through Q, and F 

smaller/greater than 1 indicates sub/super-Poissonian statis¬ 
tics. In 1 was achieved, verifying the coherent 

nature of the mechanical oscillations in their setup. For com¬ 
parison, the Poissonian statistics of a (classical) coherent state 
imply F = 1 and g^^^(O) = 1, while a thermal state would have 

g(2)(0) = 2. 

Description of the system We study the optomechanical 
setup depicted in Fig. Two optical modes a and h couple 
to a mechanical mode c via the three-mode interaction Hamil¬ 
tonian V = g{){ab^ -f a^b){c -f c^),where go is the single photon 
optomechanical coupling strength and a, b, c are the lowering 
operators of the different modes. Such an interaction has been 
implemented in Refs. EHEIEI. The optical mode h is res¬ 
onantly driven with a laser of power P , which we parametrize 
with E = yfKPjWoi) (/c is the cavity line width, and cob the 
resonance frequency of mode b). The other optical mode a is 
detuned with respect to cavity mode h and the driving laser by 
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A, and the mechanical frequency is so that the Hamilto¬ 
nian in a rotating frame for both cavities with frequency ojiy is 
H = Ho -\-V iE(b^ - b) with Hq = cOmC^c - Aa^a. 

Depending on the sign of the laser detuning, the laser either 
cools the mechanical mode (A < 0) , or gives rise to self- 
induced mechanical oscillations (A > 0). In the latter regime 
the intrinsic nonlinearity of the three-mode optomechanical 
interaction V stabilizes the mechanical oscillation at a finite 
amplitude |[3^ . We choose a detuning A = oj^ between the 
two cavities which corresponds to a resonant excitation of op¬ 
tomechanical limit cycles. In an interaction picture with re¬ 
spect to Ho the Hamiltonian is 

Hj = iE(b^ - b) goiab^c a^bc^). (2) 


We neglected here fast oscillating terms goab'^c'^ h.c., 

assuming a cavity decay rate of k ^ ojm for both cavities 
(the corrections are of order I i.e. negligible for typi¬ 
cal optomechanical crystals.). In the framework of Langevin 
equations the system dynamics is then described by 

a = -igobc^ “ f ^ (3) 

b = -igodc - ^b E ^^Kbin, (4) 

c = -igoo^b - -r ^Jycir^, (5) 


where = {bin{t)b]y)) = S(t - t') and 

{cin{t)cl^{t')) = (1 + h)d{t - t') are the two-time correlation 
functions of the Langevin noise forces. We assumed energy 
decay of the mechanical oscillator at rate y, due to coupling to 
a thermal thermal bath with mean occupation h. We adopt the 
convention from the review m that K and y are energy decay 
rates. Correspondingly, amplitudes decay at a:/2 and y/2. 

Calculation of classical amplitudes We express each of 
the operators a, b, and c as a sum of a classical (C-number) 
component and operators describing fiuctuations around it, 
such that a = a + 6a, b = (3 + 6b and c = ^ + 6c. Insert¬ 
ing this into the Langevin equations, and considering the C- 
number components only, gives rise to a coupled set of non¬ 
linear equations for the classical cavity amplitudes a and {3, 
and the (complex) mechanical amplitude f. In particular one 
finds, a = -igoPC - f and j3 = -igoc^C - ^(3+ E. We assume 
that the optical amplitudes adiabatically follow the motion of 
the mechanical oscillator which is equivalent to the conditions 
in l)y, go\(^\, go\P\ ^ K- Solving a = p = f) results in the adi¬ 
abatic solution for the optical amplitudes 




Ek 


^(f,r) = - 


Egt^C 


( 6 ) 


where = g^\C? + Inserting these optical amplitudes in 
the equation of motion for the classical mechanical amplitude 
results in ^ = -f(y -i- yoptX, where the optically mediated 
(anti)damping is 

Topt(f) — j^2 ’ 




Figure 2. (a) Optically mediated (anti)damping joptiO (bold line) as 
a function of mechanical amplitude f according to equation (|7]). The 
steady state fo is reached when 7opt(^o) = “7 (dashed line), (b) Intra 
cavity photon number \a\^ in mode a (blue dashed line), \/3\^ in mode 
b (red dash-dotted line) and total photon number 
(black solid line) is plotted as a function of mechanical ampli¬ 
tude ^ according to equation ([^. The optically induced diffusion 
Dopt = ggldcrp + \l3\^)lh^ of the mechanical oscillator scales exactly 
like the red dash-dotted line with a scale as given on the right y- 
axis. (c) Schematic phase space trajectory of the mechanical oscil¬ 
lator approaching the limit cycle attractor with amplitude ^o- In the 
co-rotating frame of the oscillator the X quadrature relates to its am¬ 
plitude and the Y quadrature to its phase. 


cf. Fig. [^. yopt is negative for all mechanical amplitudes 
and its absolute value decreases with increasing amplitude f 
according to the Lorentzian given by /z^, approaching 0 for 
f » k! go. In agreement with ll^ we define the dimensionless 
parameter 


^ lyopt(0)l 

y K^y 

which corresponds to the gain of mechanical amplification at 
zero mechanical amplitude. For < 1 the total mechanical 
damping y+yopt(0) < 0 is positive for all amplitudes, implying 
f = 0 in steady state. Above threshold, > 1, the steady 
state il = 0) is achieved for a mechanical amplitude fo such 
that yopt(fo) = cf. Fig. |^. The solution of this nonlinear 
equation is 


The solution is unique (up to the oscillator’s phase) and fully 
determined by the gain parameterside H and the single-photon 
strong-coupling parameter 2go Ik. It is instructive to contrast 
this result with the equivalent one for a conventional, two¬ 
mode (that is one mechanical and one optical mode) optome¬ 
chanical system where the mean phonon number of self in¬ 
duced limit cycles scales as the inverse of the much smaller 
ratio igolojmf instead. In view of Eq. Q it is clear that a 
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small oscillation amplitude is advantageous in order to ob¬ 
serve strong antibunching and that the three mode setup im¬ 
proves the signal approximately by a factor of 4(6 g)^/a:)^. This 
can be two orders of magnitude for typical system parame¬ 
ters of optomechanical crystals, e.g. = 217 with 

KjlTi = 500MHz and oJml'^Tr = 3.68GHz from IH. 

In the following we will set the arbitrary phase of the limit 
cycle oscillation to be zero, fo = Ifol, without loss of gener¬ 
ality. Note also that the cavity amplitudes in Eq. change 
quite significantly as the mechanical limit cycles develops, cf. 
Fig.[^, as follows from their enhanced interaction, which de¬ 
tunes the cavity from its input. 

Calculation of quantum amplitude noise The fiuctuations 
6a, 6b, and 6c with respect to these classical amplitudes fulfill 
the linearized Langevin equations 

6d = {-^6a - igo(o5b) - igoPodC + (10) 

6b = (- ^6b - igo(o6a) - igoaoSc + yfKpn, (11) 

6c = -\6c - igQ{al6b + Po6a^) + (12) 

where we consistently dropped all terms of quadratic order in 
the fiuctuations. This approximation is only valid for large 
enough amplitudes. We also introduce here the shorthand no¬ 
tation (o(o,ySo) = (^(fo) 5 y^(fo)) for the cavity amplitudes in the 
developed mechanical limit cycle. The quantum fiuctuations 
of the cavity modes can now be treated in analogy to the clas¬ 
sical amplitudes simply by setting 6a = 6h = and solving 
the resulting algebraic equation. Inserting the solutions for 6a 
and 6b back into Eq. gives the dynamics for the me¬ 

chanical mode 6c. Eor the canonical mechanical quadratures 
X = {6c + 6c^) / and F = {6c - 6c^) j cf. Fig. ||c), 
we get effective Langevin equations 

X = -\TX + VdXjv, Y = yfDYN, (13) 

with damping F, diffusion D and noise forces fulfilling 
{X^iDMO) = 6{t - f) and <7^(0, F^(t')> = - f). 

Both F = y + Fopt(^) and D = y(^ + n) + DoptiO have 
an intrinsic mechanical constant contribution and an optically 
mediated nonlinear (^-dependent) contribution. We find that 
DoptiO = at the point of the limit cycle is 

exactly as large as the vacuum contribution of the mechanical 
bath, i.e. £>opt(^o) = 5 , but Fopt(^o) = r(3 - 4/ V^) can grow 
up to three times the mechanical damping for large In total 
the damping and diffusion depicted in Fig. |^a) are at the limit 
cycle 


F(^o) = 4r(l - 1/ V^), £>(^ 0 ) = y(n + 1). (14) 

As schematically depicted in Fig. |^c), in our convention the 
T-quadrature relates to the phase of the mechanical oscilla¬ 
tor, which is subjected to undamped diffusion, cf. Eq. Gil'- 
The X-quadrature relates to the mechanical amplitude, our fo¬ 
cus of interest in this article. In particular for the phonon oc¬ 
cupation number h = c^c one finds (h) = 0(^) and 

(n^) = 0(^), such that the Fano factor is 



Figure 3. (a) Fano factor as a function of effective mechanical 

bath occupation number h and total number of photons in the cav¬ 
ity ^phk7/4go] = according to Eqs. ( [TS] ) and ( [TS] ). (b) Plot of 
fe^'HO)-l)[(goA)' ] as a function of the same parameters in units of 
the squared single-photon strong-coupling parameter (go/ k)^ accord¬ 
ing to Eqs. 03 and 03- The condition for both sub-Poissonian 
statistics (F < 1) and antibunching < 1) is visualized by the 

red contour line h - 1-2/ spR in both plots, (c) Plot of gopUfo fo^ ^P" 
timal choice of ‘F as a function of go Ik and h according to Eq. CD- 


F ^ 2{X^). Eq. gives (X^) = DjY in steady state, i.e. 
the amplitude variance is determined by the compromise of 
diffusion and effective damping, yielding for the Fano factor 


1 \+n 

2 1 - 1/V^’ 


(15) 


This is in excellent agreement with numerical results shown in 
Figure a) that were obtained by Monte-Carlo simulation of 
a master equation equivalent to the exact, nonlinear equations 
of motion in Eqs. ^ to 0- The numerics is further described 
in the Appendix 

From Eq. GD we see that for ^ » 1 the Fano factor ap¬ 
proaches (1 -f fi)/2. Therefore, we arrive at the condition 
n < 1 necessary in order to observe sub-Poissonian phonon 
statistics. For a cryogenically cooled mechanical oscillator 
h = _ 1 ) < 1 can in principle be achieved for a 

sufficiently high resonance frequency and at low temperature 
T, see |[39j|40l. However, in the present case it is possible to 
take advantage of laser cooling of the mechanical oscillator 
||4T]|42l in order to observe sub-Poissonian statistics. 

Additional Laser Cooling Consider a setup where the me¬ 
chanical oscillator is coupled to a third optical cavity of line 
width Kd which is driven below resonance such as to induce 
an additional damping of the oscillator. Eliminating this 
cooling cavity gives rise to a ‘dressed’ mechanical oscillator 
whose equation of motion is still given by ^ with an effective 
mechanical damping and occupation number 


r = ro + 7 l , 


ro^o + yini 


(16) 


ro + rL 
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Here yo is the line width and ho the occupation number of 
the bare mechanical resonance (without laser cooling), and 

= (Kdl^ojmf is the quantum limit of optomechanical laser 
cooling 14111?^ . 

In order to have F < 1 we assume laser cooling to an ef¬ 
fective phonon occupation h < 1. This comes at the cost of a 
decreased gain parameter in Eq. K = ISg^E^jK^iyo + yi), 
which can be compensated for by a somewhat more intense 
driving field. It is rather remarkable that laser cooling can help 
to observe a quantum feature such as sub-Poissonian phonon 
statistics: While laser cooling can provide a small effective oc¬ 
cupation number h ho it does so by increasing the effective 
mechanical line width y » yo by the same factor. As a result, 
the decoherence rate relevant for quantum effects, yofio = y^, 
stays constant, such that laser cooling in most cases does not 
help in order to achieve quantum effects with mechanical os¬ 
cillators. 

Experimental feasibility with current technology The re¬ 
quirements on the system parameters to have < 1 (and 

therefore E < 1) is found by inserting the mean amplitude (|^ 
and the Fano factor in the definition ([T]i of 

g(2)(0)-l =4(^f(17) 

\kI ypR-\ 

For the discussion of experimental feasibility it is more in¬ 
structive to express the gain parameter H in terms of the steady 
state total number of photons in the cavity 

«ph = \aof + \/3of = 

4go 

where we used Eqs. 0 and ([^. The circulating number of 
photons is important as it determines the heating of the me¬ 
chanical structure, which was the limiting decoherence mech¬ 
anism in recent experiments with optomechanical crystals (41 . 
In Fig.[^(a) and (b) we show the Fano factor E and g^^^(O) - 1 
(in units of g^K^) as a function of the number of photons in 
the cavity ^ph and the effective mechanical bath occupation 
number h. In view of the dependence of the Fano factor and 
the second order coherence function on cf. Eqs. ( p~5] ) and 
(rf\ respectively, it is clear that there is an optimal number 
of circulating photons minimizing g^^^(O) for given h and sin¬ 
gle photon strong coupling parameter go/ k. The minimum is 
reached at n^^ [A:y/4gQj = (3 -r h)/(l - h) and is given by 

which is illustrated in Fig. [^. Thus, a large single-photon 
coupling helps, but is not strictly required, to create a ro¬ 
bust signal to verify antibunching. We conclude that a sub- 
Poissonian phonon laser state can be prepared and verified 
outside the single-photon strong-coupling regime and for 
small but finite effective (cf. Eq. ( p^ ) bath occupation h 
by detecting photon antibunching in the reflected light. We 
emphasize that phonon antibunching can be observed already 
in a regime of few circulating photons ^ph ^ 1 • 



Figure 4. (a) Comparison of analytical results from Eq. GD to nu¬ 
merical results for Fano factor F with increasing go/K = 0.25,0.5,1. 
(square, diamond, circle). The parameter goE/K^ = 0.04 is fixed to 
stay well inside the regime of validity of the adiabatic elimination. In 
this plot n = 0 but for finite temperature the agreement of numerics 
with Eq. GD is equally good, (b) Negativity of the Wigner func¬ 
tion W of the mechanical oscillator in steady state calculated with 
QuTiP’s steady state solver. We define the negativity as the quotient 
of the smallest and the largest value of W. In this plot the bath oc¬ 
cupation is w = 0.25,0.5,1,2 for the increasing curves. The driving 
field E - 0.07/c is constant in each of these plots, to stay well in the 
regime of 1 for numerical simplicity. Each point is optimized 

over K by varying y. 


Readout The readout of the - possibly antibunched - 
phonon statistics can be implemented in analogy to CD using 
the cooling laser mode d. In the sideband resolved (/cj cOm) 
and linear (g^f co^) regime the dynamics of laser cooling 
can be understood as a continuous coherent state swap inter¬ 
action cdf c^d Eia. The phonon statistics of d can then be 
measured by counting the photons in the output of the cooling 
cavity d at the sideband frequency +0Jm CD. Hence with this 
readout scheme phonon antibunching is detected via photon 
antibunching. 

Experimental case study Currently the highest reported 
value for the coupling in optomechanical crystals is go/^/r = 
1.1 MHz (4311 . The lowest cavity decay rate in a photonic 
crystal is, to our knowledge, /c = 20 MHz (44l . While the 
best ratio achieved in a single device is go/K = 0.007 0, 
combining the best values in one device would already reach 
go/K ^ 0.055. The lowest reported effective bath occupation 
reached with optomechanical cooling is = 0.85 0, using 
a dilution refrigerator mechanical oscillators have even been 
cooled down below h < 0.07. Assuming a slightly more op¬ 
timistic go/K = 0.1, an effective environmental temperature 
of 200 mK and a mechanical frequency of 5GHz the devia¬ 
tion of gopt(^) 1 according to Eq. ( p^ will be 2.5 per 
mille. Further improvements on go and k are expected using 
new designs and fabrication methods, so that reaching a sig¬ 
nal of gopt(O) - 1 on the order of a few per cent is a realistic 
prospect for the near future, cf. Fig. 

Outlook: Towards the single-photon strong-coupling 

regime Our linearized model is strictly valid only for go /k 
1. We can however expect qualitative agreement to some ex¬ 
tend even for larger go/K. The deviation of E from equation 
GU in this regime are plotted in Fig. Ef- Strongly sub- 
Poissonian states with small limit cycle amplitude {h) fea¬ 
ture a negative Wigner function (23l. As discussed above 
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{h) ~ (a:/^o)^ • It is therefore reasonable to expect negative 
mechanical Wigner density with g^jK approaching the single¬ 
photon strong-coupling regime. Indeed we numerically find 
that negative Wigner density is possible for larger g^jK as de¬ 
picted in Fig. 1^. All numerical calculations were done with 
QuTiP HSlllSl, the details of the methods are discussed in the 
Appendix. 

Conclusion Using an optomechanical setup with two op¬ 
tical modes brings experimental demonstration of both sub- 
Poissonian phonon statistics and optomechanically induced 
phonon and photon antibunching in reach of today’s technol¬ 
ogy. For parameters approaching the single-photon strong¬ 
coupling regime the limit cycle states can even feature a neg¬ 
ative mechanical Wigner function. 
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numerics The steady state of the system was calculated 
using QuTiP 14511461 . For Fig. |^, where the mechanical 
amplitudes are small due to the large go the Hilbert space 
has moderate size and we used a direct steady state solver for 
density matrices. For Fig. the Hilbert space is (in gen¬ 
eral) too large for this and we had to use Monte-Carlo trajec¬ 
tories ll47Vl49ll for the wave function and average over many 
runs to obtain a density matrix. Each trajectory 10^/0) had a 
coherent state with random, independent and identically dis¬ 
tributed Gaussian amplitudes, ~ 1) around the ana¬ 


lytical steady state amplitude f from Eq. as initial state 
for the oscillator and vacuum as initial state for both optical 
modes. The system was then evolved for a time t = 5ly with 
Hamiltonian Q and the Lindblad operator L = Lc-^-L^, where 
Lcp = Kapa^ap-^pa^a+Kbpb^-^b^bp-^pb^b dindLm = 
ycpc^ -^c^cp- ^pc^c. The calculation was done in a displaced 
frame around the mean amplitude of the mechanical oscillator 
and cavity modes. We then used cr = \i//j(T)\ to cal¬ 
culate {h}cr and which in turn gives the Fano factor F 

and g^^^(O). The time r was chosen such that both mean val¬ 
ues had already relaxed to steady state, compare the damping 
in equation ( p3] ), while the phase has still not diffused away 
too far from f so that a small Hilbert space around the mean 
mechanical amplitude was sufficient for the simulation. 


